High— Order Unstructured Lagrangian 

One— Step WENO Finite Volume Schemes for 

Non— conservative Hyperbolic Systems: 

Applications to Compressible Multi— Phase 

Flows 

en 

o 

(N 

5_i Michael Dumbser ^ Walter Boscheri ^ 

Oh 
<^ ^Laboratory of Applied Mathematics 

Department of Civil, Environmental and Mechanical Engineering 

University of Trento, Via Mesiano 11, 1-38123 Trento, Italy 



l> 



< 

^ Abstract 



cd 



In this article we present the first better than second order accurate unstructured 

Lagrangian-type one-step WENO finite volume scheme for the solution of hyper- 

I I bolic partial differential equations with non-conservative products. The method 

achieves high order of accuracy in space together with essentially non-oscillatory 

<^ behaviour using a nonlinear WENO reconstruction operator on unstructured tri- 

^<D angular meshes. High order accuracy in time is obtained via a local Lagrangian 

'"^ space-time Galerkin predictor method that evolves the spatial reconstruction poly- 

■^ nomials in time within each element. The final one-step finite volume scheme is 

.^ derived by integration over a moving space-time control volume, where the non- 

^^ conservative products are treated by a path-conservative approach that defines 

^^ the jump terms on the element boundaries. The entire method is formulated as 

. . an Arbitrary-Lagrangian-Eulerian (ALE) method, where the mesh velocity can be 

. ^ chosen independently of the fluid velocity. 

rN The new scheme is applied to the full seven-equation Baer-Nunziato model of 

C^ compressible multi-phase flows in two space dimensions. The use of a Lagrangian 

approach allows an excellent resolution of the solid contact and the resolution of 
jumps in the volume fraction. The high order of accuracy of the scheme in space and 
time is confirmed via a numerical convergence study. Finally, the proposed method 
is also applied to a reduced version of the compressible Baer-Nunziato model for 
the simulation of free surface water waves in moving domains. In particular, the 
phenomenon of sloshing is studied in a moving water tank and comparisons with 
experimental data are provided. 

Key words: Arbitrary-Lagrangian-Eulerian (ALE) scheme, WENO finite volume 
scheme, path-conservative scheme, unstructured meshes, high order in space and 
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1 Introduction 



Multi-phase flow problems, such as liquid-vapour and solid-gas flows are en- 
countered in numerous natural processes, such as avalanches, meteorological 
flows with cloud formation, volcano explosions, sediment transport in rivers 
and on the coast, granular flows in landslides, etc., as well as in many industrial 
applications, e.g., in aerospace engineering, automotive industry, petroleum 
and chemical process engineering, nuclear reactor safety, paper and food man- 
ufacturing and renewable energy production. Most of the industrial appli- 
cations are concerned with compressible multi-phase flows as they appear 
for example in combustion processes of liquid and solid fuels in car, aircraft 
and rocket engines, but also in solid bio-mass combustion processes. Already 
the mathematical description of such flows is quite complex and up to now 
there is no universally agreed model for such flows. One wide-spread model 
is the Baer-Nunziato model for compressible two-phase flow, which has been 
introduced by Baer and Nunziato in [7] for detonation waves in solid-gas 
combustion processes and which has been further extended by Saurel and 
Abgrall to liquid-gas flows in [119]. It is therefore also often called in lit- 
erature the Saurel-Abgrall model. The main difference between the original 
Baer-Nunziato model and the Saurel-Abgrall model is the deflnition of the 
pressure and velocity at the interface. In the present paper, we will use the orig- 
inal choice of Baer-Nunziato, which has also been used in several papers about 
the exact solution of the Riemann-Problem of the Baer-Nunziato model, see 
[5,39,120]. A reduced flve-equation model has been proposed in [102] and 
approximate Riemann solvers of Baer-Nunziato-type models of compressible 
multi-phase flows can be found for example in [51,55,129,127]. 

The full seven-equation Baer-Nunziato model with inter-phase drag and pres- 
sure relaxation is given by the following non-conservative system of partial 
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differential equations: 
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(0lPl) + V ■ (0lPlUi) = 0, 

(0lPlUi) + V ■ (0iPiUiUi) + V0iPi = P/V01 - A (Ui - U2) , 

{(piPiEi) + V ■ {{(pipiEi + (pipi) ui) = -pidtcj)! - Aui ■ (ui - U2) , 

:P2) + V ■ (02P2U2) = 0, ^ (1) 

;P2U2) + V ■ (02P2U2U2) + V02P2 = P/V02 - A (U2 - Ui) , 
:p2E2) + V ■ ((02P2^2 + 02^2) U2) = J5/(9t0i - A Uj ■ (U2 - Ui) , 
UiV0l = U{pi -P2)- 



In the entire article, the system is closed by the so-called stiffened gas equation 
of state (EOS) for each phase: 

Pk + IkTTk .„^ 

efc = —. 7T- 2 

Pkilk - 1) 

Here, (pk denotes the volume fraction of phase k, pk is the density, Uk is the 
velocity vector, E^ = ek + |uk^ and e^ are the phase specific total and internal 
energies, respectively, A is a parameter characterizing the friction between both 
phases and u characterizes pressure relaxation. For consistency the sum of the 
volume fractions must always be unity, i.e. 0i + 02 = 1- In the literature, 
one of the phases is often also called the solid phase and the other one the 
gas phase. Defining arbitrarily the first phase as the solid phase in the rest of 
the paper we will therefore use the subscripts 1 and s as well as 2 and g as 
synonyms. For the interface velocity and pressure ui and pi we choose ui = ui 
and Pi = p2 respectively, according to [7], although other choices are possible, 
see e.g. the paper by Saurel and Abgrall [119]. We can cast system (1) in the 
general non-conservative form (3) below 

^ + V-F(Q) + B(Q)-VQ = S(Q), xGficR^tGR^ (3) 

where Q = (0ipi, 0iPiUi, (pipiEi, 02P2, 02P2U2, (P2P2E2, 0i) e Qq C W is 
the state vector, F = (f , g) is the flux tensor, i.e. the purely conservative part 
of the PDF system, B = (Bi, B2) contains the purely non-conservative part of 
the system in block-matrix notation and S(Q) is the vector of algebraic source 
terms, which, in our case, contains the inter-phase drag and the pressure 
relaxation terms. 

We furthermore introduce the abbreviation P = P(Q, VQ) = B(Q) ■ VQ. 
System (3) can also be written equivalently in the following quasi-linear form 

^ + A(Q)-VQ = S(Q), (4) 



with A(Q) = (Ai, A2) = 9F(Q)/9Q + B(Q). The fluxes F, the system ma- 
trix A(Q) and the source term vector S(Q) can be readily computed from (1). 
Hyperbolicity of the Baer-Nunziato model and exact Riemann solvers have 
been studied in [5,39,120]. A very important requirement for numerical meth- 
ods used to solve the Baer-Nunziato model is the sharp resolution of material 
interfaces, i.e. jumps in the volume fractions (pk- Improved resolution of the 
material interfaces can be achieved using high order finite volume methods 
together with little diffusive Riemann solvers, such as the Osher scheme and 
the HLLC method, both of which have already been successfully applied to 
the Baer-Nunziato model in [55] and [129], respectively. However, if the inter- 
face velocity is small compared to the sound speed of at least one of the two 
phases, significant numerical smearing of the material contact will occur even 
with high order schemes and little diffusive Riemann solvers due to the small 
CFL number associated with the material wave. It may therefore be necessary 
to further improve the resolution of material interfaces by using a Lagrangian 
method instead of an Eulerian one, since in the Lagrangian case the mesh 
moves with the flow field and therefore allows a sharp tracking of material in- 
terfaces independent of the CFL number associated with the material contact 
wave. 

The significantly improved resolution of material interfaces led to intensive 
research on Lagrangian schemes in the past decades. The construction of La- 
grangian schemes can start either directly from the conservative quantities 
such as mass, momentum and total energy [95,124], or from the nonconserva- 
tive form of the governing equations, as proposed in [14,20,136]. Furthermore 
existing Lagrangian schemes can be divided into staggered mesh and cell- 
centered approaches, or combinations of them [90]. Cell centered Godunov- 
type finite volume schemes together with the first Roe linearization for La- 
grangian gas dynamics have been proposed by Munz in [101]. Munz found that 
in the Lagrangian framework, the Roe-averaged velocity for the equations of 
gas dynamics is simply given by the arithmetic mean, while in the Eulerian 
frame the Roe average is a more complicated function of the left and right 
densities and velocities. A cell-centered Godunov scheme has been proposed 
by Carre et al. [21] for Lagrangian gas dynamics on general multi-dimensional 
unstructured meshes and in [40] Despres and Mazeran introduce a new formu- 
lation of the multidimensional Euler equations in Lagrangian coordinates as 
a system of conservation laws associated with constraints. Furthermore they 
propose a way to evolve in a coupled manner both the physical and the ge- 
ometrical part of the system [41], writing the two-dimensional equations of 
gas dynamics in Lagrangian coordinates together with the evolution of the 
geometry as a weakly hyperbolic system of conservation laws. This allows the 
authors to design a finite volume scheme for the discretization of Lagrangian 
gas dynamics on moving meshes, based on the symmetrization of the formu- 
lation of the physical part. In a recent work Despres et al. [36] propose a 
new method designed for cell-centered Lagrangian schemes, which is transla- 



tion invariant and suitable for curved meshes. General polygonal grids have 
been considered by Maire et al. [94,92,91], who develop a general formalism 
to derive first and second order cell-centered Lagrangian schemes in multiple 
space dimensions. By the use of a node-centered solver [94], the authors obtain 
the time derivatives of the fluxes and hence a second order method in space 
and time. Lagrangian schemes for multi-material flows have been successfully 
introduced in [18,15,19] and Lagrangian schemes with additional symmetry 
preserving properties in cylindrical geometries have been proposed for exam- 
ple in [33,34,93,89]. A multi-scale cell centered Godunov-type finite volume 
scheme for Lagrangian hydrodynamics has been introduced in [96]. All the 
Lagrangian schemes listed before are at most second order accurate in space 
and time. 

Cheng and Shu were the first to introduce a third order accurate essentially 
non-oscillatory (ENO) reconstruction operator into Godunov-type Lagrangian 
finite volume schemes [32,87]. Their cell centered Lagrangian finite volume 
methods achieve also high order in time, either by using the method of lines 
(MOL) approach based on a third order TVD Runge-Kutta time discretiza- 
tion, or by using a high order one-step Lax-Wendroff-type time stepping. 
Higher order unstructured Lagrangian finite element methods have been re- 
cently investigated in [108,121]. In a very recent paper Dumbser et al. [56] 
propose a new class of high order accurate Lagrangian-type one-step WENO 
finite volume schemes for the solution of stiff hyperbolic balance laws. In [16] 
this class of high order one-step Lagrangian WENO finite volume schemes has 
been extended to unstructured triangular meshes for the conservative case and 
for single-phase fiows. The present paper is concerned with its extension to 
non-conservative systems and applications to compressible multi-phase fiows. 
To our knowledge, the work presented in this article is the first better than 
second order accurate unstructured Lagrangian-type finite volume scheme for 
non-conservative hyperbolic systems with applications to compressible multi- 
phase fiows. 

To round-off this introduction, we briefiy refer also to other Lagrangian- 
type schemes and alternative Eulerian schemes with improved resolution of 
material interfaces. The following list of references does not pretend to be 
complete. Among alternative Lagrangian schemes one should also mention 
for example meshless particle schemes, such as the smooth particle hydro- 
dynamics (SPH) method [98,64,65,66], which has been successfully used to 
simulate fiows in complex domains with large deformations. Since the SPH 
approach is a fully Lagrangian method the mesh moves with the local fiuid 
velocity, whereas in Arbitrary Lagrangian Eulerian (ALE) schemes, see e.g. 
[74,113,124,42,63,62,27], the mesh moves with an arbitrary mesh velocity that 
does not necessarily coincide with the real fiuid velocity. This adds an extra 
degree of fiexibility and generality to the scheme. The classical SPH method 
is a truly meshless scheme, but exhibits also a series of problems, such as the 



need for artificial stabilization terms and the lack of zerotli order consistency. 
These problems have been overcome by the so-called Particle Finite Element 
Method (PFEM), which is a Lagrangian (or Arbitrary-Lagrangian-Eulerian) 
finite element scheme on moving point clouds, [78,114,107,84,79,106], and has 
been successfully applied to incompressible multi-material flows and fluid- 
solid interaction problems. Furthermore, the reader will also find so-called 
Semi-Lagrangian schemes in literature, which are mainly used for solving 
transport equations [118,71]. Here, the numerical solution at the new time 
is computed from the known solution at the current time by following back- 
ward in time the Lagrangian trajectories of the fluid to the foot-point of the 
trajectory. Semi-Lagrangian schemes are therefore nothing else than a numer- 
ical scheme based on the method of characteristics. Since in general the end- 
point does not coincide with a grid point, an interpolation formula is required 
in order to evaluate the unknown solution, see e.g. [25,26,86,77,115,38,17]. 
Note that in Semi-Lagrangian algorithms the mesh is fixed, like in a classi- 
cal Eulerian methods. An alternative to Lagrangian methods for the accurate 
resolution of material interfaces has been developed in the Eulerian frame- 
work on fixed meshes under the form of ghost-fluid and level-set methods 
[60,61,67,109,100], as well as the volume of fluid method [75,117,88]. 

The rest of the paper is organized as follows: in Section 2 the high-order 
path-conservative Arbitrary-Lagrangian-Eulerian (ALE) one-step WENO fi- 
nite volume scheme is described on unstructured triangular meshes and a 
numerical convergence study on a smooth unsteady test problem is carried 
out, to assess the designed order of accuracy of the scheme in space and time. 
In Section 3 some classical academic benchmark problems with exact or quasi- 
exact reference solution are solved to verify the robustness of the method in 
the presence of shock waves and significant mesh distorsion. Finally, in Sec- 
tion 4 our new high order path-conservative ALE method is applied to a 
free-surface flow problem concerning sloshing in a moving tank. For this pur- 
pose, a reduced version of the Baer-Nunziato model is solved, which has been 
introduced for the simulation of free surface flows on fixed grids in [46,44] . The 
paper is rounded-off by some concluding remarks and an outlook to future 
research in Section 5. 



2 Numerical Method 



The time-dependent computational domain is denoted by Q{t) C R^ and is 
discretized at a given time t" by a set of conforming triangles T", the union 
of which is the current triangulation T^ of the domain Vtit^) = 1]" and can 



be expressed as 

Ne 

r^ = [j Tp. (5) 

i=l 

Here, Ne is the number of elements used to discretize the domain. 

On triangular meshes it is convenient to introduce also a local spatial reference 
coordinate system C, — tj, which maps the physical element T" in the current 
configuration to the reference element denoted by Tg. The spatial mapping 
reads 

X = x(^, n = XI, + (X^,, - XI,) e + (X^,, - XI,) V (6) 

where $, = {C,, rj) and x = (a;, y) are the vectors of the spatial coordinates in the 
reference system and the physical system, respectively, and X^j = [XJl^^Y^^ 
is the vector of physical coordinates of the /c-th vertex of triangle T" at time 
t". The unit triangle Tg is defined by the nodes ^g i = {.ie^i^Ve,!) = (0,0), 

L,2 = fe,2,^e,2) = (1,0) and ^g_3 = (Ce,2,^e,2) = (0,1). 

The cell averages, which represent the data that are stored and evolved in 
time within a finite volume scheme, are defined at time t" as usual by 



Qr = ^^/ Q(x,r)rfx. (7) 

Here, |T"| denotes the area of triangle T". Higher order in space can be 
achieved by reconstructing piecewise higher order polynomials w/i(x, t") from 
the cell averages defined above in Eqn. (7). For this purpose, we employ 
a higher order WENO reconstruction procedure following [82,52,53]. Other 
WENO schemes on structured and unstructured meshes can be found, e.g., in 
[9,80] and [68,128,133,76,138,3], respectively. Instead of a WENO method, al- 
ternative high order nonlinear reconstruction operators can be used as well, see 
e.g. [125,1,35]. For the convenience of the reader, the component-wise WENO 
reconstruction procedure is briefly summarized in the next section, for details 
we refer to [52] . A possible polynomial WENO reconstruction in characteristic 
variables can be found in [53]. 



2.1 WENO Reconstruction 



The reconstructed solution w/j(x, t") is represented by piecewise polynomials 
of degree M and is obtained from the given cell averages in an appropriate 
neighborhood of element T", called the reconstruction stencil 5/. The number 
of elements inside each reconstruction stencil is denoted by Ue- In 2D we use 
in total 7 reconstruction stencils for each element, one central stencil, three 
forward stencils and three backward stencils, as proposed in [82,52], hence 
1 < s < 7. According to [12] the total number of stencil elements Ue must be 



larger than the number of degrees of freedom Ai = [M + 1)(M + 2)/2 of a 
polynomial of degree M. Typically we take rig = 2A^ in two space dimensions. 

The reconstruction polynomial for each candidate stencil s for triangle T" is 
written in terms of spatial basis functions ipii^,) as 



w^(x,r) = E^KI)wi:r:=^z(|)w"'^ 



1=1 

where the mapping x = x(^, t"') is given by Eqn. (6). In the rest of the paper 
we will use classical tensor index notation based on the Einstein summation 
convention, which implies summation over two equal indices. The number 
of the unknown degrees of freedom to be reconstructed for each element is 
Ai. The basis functions ipi{^) are the orthogonal basis functions described in 
[43,81,37]. 

The reconstruction on each stencil S^ is based on integral conservation, i.e. 

^ / M^)^?fd^ = Q ■ , vt; g St. (9) 

j 

Since n^ > Ai the above system (9) is an over-determined linear algebraic 
system that is solved for w"^* using a constrained least-squares technique, see 
[52]. The linear constraint is that Eqn. (9) holds exactly at least for element 
T^. The multi-dimensional integrals appearing in the expression above are 
evaluated using Gaussian quadrature formulae of suitable order, see [126] for 
details. Since the triangles are moving in a Lagrangian scheme, the small linear 
systems (9) are solved for each element at the beginning of each time step. 
However, the choice of the stencils S^ remains fixed for all times. 

To obtain a higher order essentially non-oscillatory polynomial the scheme 
must be nonlinear, in order to circumvent the Godunov theorem that states 
that linear monotone schemes are at most of order one. The final nonlinear 
WENO reconstruction polynomial is therefore computed in the following way. 
First, the smoothness of each reconstruction polynomial obtained on stencil 
S^ is measured by a so-called oscillation indicator cr^ [80]. According to [52] 
the smoothness indicator can be easily computed on the reference element 
using the (universal) oscillation indicator matrix S^^ as follows: 

<Ts = ^Im^l, W^„ S,™ = )_^ / d^dr]. (10) 



a+l3<M rp 

The nonlinear weights ujg are defined by 



d^aQ^p d^aQ^p 



(o-s + e) Y.a^q 



where we use e = 10~^^, r = 8, As = 1 for the one-sided stencils and A^ = 10^ 
for the central stencil, according to [52]. The final nonlinear WENO recon- 
struction polynomial and its coefficients are then given by 

M 

w,(x,r)=Ev^K^K„ with w[:, = ^u.,w[:r. (12) 

/=1 s 



In order to reduce the computational cost associated with the nonlinear WENO 
reconstruction procedure outlined above, a high order one-step time dis- 
cretization is used so that reconstruction has to be performed only once for 
each time step. 



2.2 Local Space-Time Predictor on Moving Curved Meshes 



Higher order of accuracy in time is achieved by an element-local predictor 
stage that evolves the reconstructed polynomials w/i(x, t") locally in time 
within each element Tj(t) during the time interval [t"; t""*"^], see [50,47,51,72,69]. 
Such an element-local time-evolution procedure has also been used within the 
MUSCL scheme of van Leer [134] and the original ENO scheme of Harten et 
al. [70], who called this element-local predictor with initial data Wft(x, t") the 
solution of a Cauchy problem in the small, since no information from neighbor 
elements is used. The coupling with the neighbor elements occurs only later 
in the final one-step finite volume scheme. 

The local data evolution step leads for each element to piecewise space-time 
polynomials of degree M, denoted by q/i(x, t) in the following. While the 
original ENO scheme of Harten et al. uses a higher order Taylor series in 
time together with the strong differential form of the PDE to substitute time- 
derivatives with space derivatives (the so-called Cauchy-Kovalewski or Lax- 
Wendroff procedure [85]), here a weak formulation of the PDE in space-time 
is used. The resulting method does not require the computation of higher 
order derivatives, but just pointwise evaluations of the fluxes, source terms 
and non-conservative products appearing in the PDE. Such a local space- 
time Galerkin predictor scheme was introduced for the Eulerian framework 
in [50,47,51,72] and is extended here to the Lagrangian framework on moving 
curved space-time elements for PDE with nonconservative products. 

Let X = (x, y, t) denote the physical space-time coordinate vector and ^ = 
(C) 77, r) the reference space-time coordinate vector, while x = (x, y) and ^ = 
(^, rj) are the purely spatial coordinate vectors already introduced previously. 
Let furthermore 61 = Oi{$,) = 9i{C,,ri,T) be a space-time basis function deflned 
by the Lagrange interpolation polynomials passing through the space-time 
nodes $,^ = {^myVmjTm) specifled according to [47] and also depicted for the 




Fig. 1. Iso-parametric mapping of the space-time reference element (left) to the 
physical space-time element (right) used within the local space-time Galerkin pre- 
dictor. 

case M = 2 in Figure 1. The Lagrange interpolation polynomials define a 
nodal basis which satisfies the interpolation property 



Oliir. 



Hm^ 



(13) 



with the usual Kronecker symbol 5im- The element-local space-time predic- 
tor solution q/i, the fluxes F/i = (fft,,g/i), the source term S/i and the non- 
conservative product Pfe = B(q/i) ■ Vq/, are approximated within the space- 
time element Ti{t) x [t";t'^+i] as 



Sh = S,(|) = ei{i) S;,, 



p, = p,,(|) = ^KI)Pm- (14) 



Since a nodal basis is used, the degrees of freedom for F/j, Sh and P/i can be 
simply computed pointwise from q/j as 

Fi^i = F{qi^i), S/,i = S(q/,i), Pi^i = P(q/,i, Vq^^j), Vq^ ^ = V6'^(|i)qm,i. 

_ - (15) 

The degrees of freedom Vq^ j represent the gradient of q^ in node ^;. 

In the present paper an isoparametric mapping is used between the physical 
space-time coordinate vector x and the reference space-time coordinate vector 
^, i.e. the mapping is also represented by the nodal basis functions 6i. Hence, 



^{^) = 9i{^)5li, 



m = 9i{^)ti, 



(16) 



where the degrees of freedom x;,j = {xi^i,yi^i) denote the in general unknown 
vector of physical coordinates in space of the moving space-time control vol- 
ume and the ti denote the known degrees of freedom of the physical time at 



10 



each space-time node x^^j = {xi^i,yi^i,ti). The general representation of the 
mapping in time given above in Eqn. (16) simphfies to 



6 ir 



rAt, 



r 



t-r- 

At 



L] Lr. 



nAt, 



(17) 



where t" is the current time and At is the time step. Hence, t^ = t^, = and 

tr = At. The isoparametric mapping (16) allow us to transform the physical 

space-time element to the unit reference space-time element Tg x [0,1]. A 

sketch of this mapping is depicted in Figure 1 and the Jacobian matrix of the 

transformation reads 

/ \ 



J, 



St 



dx 

di 



y( Vv Vt 
At 



Its inverse is given by 



/-I 

'^st 



di 
dx 


1 \ 

ix iy it 

Vx Vy Vt 




[o o^J 



(19) 



The local reference system and the inverse of the associated Jacobian matrix 
(19) are used to rewrite the governing PDE (3) as 



dQ 
using r^ 



At 



i-i-(ir-^-B(Q,.(i)%.Q- 



AtS(Q), 



and Tt = ^ according to (17) and with the notation 




dx 

A. 
dy 



ix Vx 
Zy Vy 



(a 



9x 



V<. 



(20) 



(21) 



The term -^ ■ ^ is due to the Lagrangian mesh motion and is zero in the 
Eulerian case, i.e. for fixed meshes. 



We further introduce the abbreviation 

T 



H = AtS - At 



9Q di 

di ' dt 



di 

(9x 



V^-F + B(Q) 



and its numerical approximation by 

H/i = 6i{$,)'iii^i, 



'de 

(9x 



V^Q 



(22) 



(23) 
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as well as the two operators 

1 
[/, qY = Jf(^^ ^)9it r)d^dv, (/, g) = JJf{t T)g{t T)d^dvdT (24) 

Te Te 

that denote the scalar products of two functions / and g over the spatial 
reference element Tg at time r and over the space-time reference element 
Tg X [0, 1], respectively. 

Multiplication of (20) with space-time test functions 0k{$,), integration over 
the space-time reference element Te x [0, 1] and inserting (14) and (23) yields 

The term on the left hand side can be integrated by parts in time, which also 
allows to introduce the initial condition of the local Cauchy problem in a weak 
form as follows: 

Mi, 1), di{i, l)Y qz, - /^, 6\ q,, = [4(^, Q),Ui)f ^h + {Ou, 6i) %,. 

(25) 
With the definitions 

(26) 
the above expression, which is a nonlinear algebraic equation system for the 
unknown coefficients q^^j can be written in a more compact matrix form as 

i^iq,,, = Fowj;, + MH,,„ (27) 

and is conveniently solved using the following iterative scheme: 

q[f = K^^ (Fow[;, + MH[j , (28) 

where r denotes the iteration number. For an efficient initial guess based on 
a second order MUSCL-type scheme, see [72]. 

Since the mesh is moving we also have to consider the evolution of the vertex 
coordinates of the local space-time element, whose motion is described by the 
ODE system 

^ = V(Q,x,t), (29) 

where V = V(Q,x, t) is the local mesh velocity. In the Arbitrary Lagrangian- 
Eulerian (ALE) framework used in this article, the mesh velocity can be chosen 
independently of the fluid velocity, hence for V = the scheme reduces to 
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a pure Eulerian approach, but if V coincides with the local fluid velocity 
V = v(Q) one obtains a Lagrangian method. The discrete velocity field inside 
element Tj(t) can be expressed as 



yh = 9iitr)yi„ (30) 



with Yi^i = V{qi^i,xi^i,ii). 



As in [56] the ODE (29) can be solved for the unknown coordinate vector x; 
by using again the local space-time DG method: 

i^iXM = [4(^,0), x(^,r)]° + AtMVz,„ (31) 

with x(^, t") given according to the spatial mapping (6) based on the known 
vertex coordinates of triangle T" at time t". This results in the following 
iteration scheme for the element-local space-time predictor for the nodal co- 
ordinates: 

x[,f = K^' {MtO)Mtn]' + ^tMVl) . (32) 

Eqn. (32) is iterated together with Eqn. (28). The iteration stops when the 
residuals of (28) and (32) are less than a prescribed tolerance, which we set 
to 10^^^ for all examples shown below. 

Once we have carried out the above procedure for all the elements of the 
computational domain, we end up with an element-local predictor for the 
numerical solution q/j, as well as for the mesh velocity V/j. 

Next, we have to update the mesh globally. Let us denote with Vk the neigh- 
borhood of vertex number k, i.e. all those elements that have in common the 
node number k. The number of elements in the neighborhood Vk is denoted 
with Nk- Since the velocity of each vertex is defined by the local predictor 
within each element, one has to deal with several, in general different, veloci- 
ties for the same node, since all elements belonging to Vk will in general give 
a different velocity contribution, according to their element-local predictor. 
Since we do not admit the geometry to be discontinuous, we decide to fix a 
unique node velocity Vj^ to move the node. The final velocity is chosen to 
be the average velocity considering all the contributions V^ ■ of the vertex 
neighborhood as 

^l = W ^ ^l^^ ^^^^ ^.^- = i hi(^eMk)^^)dA V,,,. (33) 

The ^e,m(fc) s-re the vertex coordinates of the reference triangle Tg correspond- 
ing to vertex number k, hence m = m{k) with 1 < tti < 3 is a mapping from 
the global node number k to the element-local vertex number. Since each 
node now has its own unique velocity, the vertex coordinates can be moved 
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according to 

Xri=X^ + AtV:, (34) 

and we can update all the other geometric quantities needed for the compu- 
tation, e.g. normal vectors, volumes, side lengths, barycenter positions, etc. 

2.3 Path- Conservative One-Step Finite-Volume Scheme 

The governing PDE system (3) can be rewritten more compactly using the 
following space-time divergence form 

V-F + B(Q)-VQ = S(Q), (35) 

with the space-time nabla operator 

and the space-time flux tensor and system matrices 

F = (f, g, Q), B = (Bi,B2,0), A = |^ + B. (37) 

With (36) and (37) the quasi-linear form of the PDE (35) reads 

A(Q)-VQ = S(Q). (38) 

Integration over a space-time control volume C" = Ti(t) x [t"';t"+-'^] yields 

/" V ■ F d^dt + /" B(Q) ■ VQ d^dt = /" S(Q) dxdt. (39) 

Application of the theorem of Gauss allows us to write the first space-time 
volume integral on the left as a flux integral over the space-time surface (9C". 
Furthermore, the non-conservative product is integrated by using a path- 
conservative approach [132,112,111,23,99,22,116,49,51,55], which follows the 
theory of Dal Maso-Le Floch and Murat [97] and defines the non-conservative 
term as a Borel measure. For the known limitations and deficiencies of path- 
conservative schemes see [24,2]. Thus 

/ (f + d) ■ firfS + f B(Q) ■ VQ d^dt = I S(Q) rfxrft, (40) 

where n = {fix, Uy, fit) is the outward pointing space-time unit normal vector 
on the space-time surface dC^, and D is a term that takes into account 
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potential jumps of Q on the element boundaries according to the path integral 



1 



t ~ / \ (7W 

D ■ fi = y B (*(Q-, Q+, s)) ■ fi ^ ds. (41) 




Throughout the entire paper and according to [111,23,51,55] we use the fol- 
lowing straight-line segment path 

* = *(Q-, Q+, s) = q, + s(Q+ - Q-), (42) 

for which the jump term above (41) simplifies to 



(43) 



The space-time surface dC^ above involves overall five space-time sub-surfaces, 
as depicted in Figure 2: 

9Cr = I U 9^5 I U T;- U T:^\ (44) 

where Mi denotes the so-called Neumann neighborhood of triangle Tj(t), i.e. 
the set of directly adjacent triangles Tj(t) that share a common edge dTij(t) 
with triangle Tj(t). The common space-time edge (9C"- during the time interval 
[t";t"+^] is denoted above by dC^^ = dTij{t) x [t";r+i]. 

The upper space-time sub-surface T""*"^ and the lower space-time sub-surface 
T" are parametrized byO<.^<lAO<r7<l— ,^ and the mapping (6). They 
are orthogonal to the time coordinate, hence for these faces the space-time 
unit normal vectors simply read n = (0,0,1) for T"^^ and fi = (0,0,-1) 
for T", respectively. The lateral space-time sub-faces dC^ are defined using a 
simple bilinear parametrization, since the old vertex coordinates X"^ are given 
and the new ones X"^"*"^ are known from (34). 



ac;;. = i(x,r) = ^/3,(x,r)X" , 0<x<l, 0<r<l, (45) 



where (x, t) represents a side-aligned local reference system according to Fig- 
ure 2. The X";, are the physical space-time coordinate vectors for the four 
vertices that define the lateral space-time sub-surface dC^y If X^-^^ and X^-g 
denote the two spatial nodes at time t" that define the common spatial edge 
dTijit"'), then the four vectors X^-;, are given by 
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(46) 



The /3fc(x, t) are a set of bilinear basis functions, which are defined as 

/3i(x, r) = (1 - x)(l - r), (S^ix, r) = x(l - r), 



/33(X, r) = XT, f3,{x, r) = (1 - x)r. (47) 

From (46) and (47) it follows that the temporal mapping is again simply 



t = t" + r At, hence ty = and t^ = At. The determinant of the coordinate 







T^7l - 




(a) 



(b) 



Fig. 2. Physical space-time element (a) and parametrization of the lateral space — 
time subsurface dCf", (b). 

transformation and the resulting space-time unit normal vector fijj of the 
sub-surface dC^j can be computed as follows: 



\dcrA 



dx dr 



n 



V 



|^f)/l^-;i- 



(4J 



Note that the integration over a closed space-time volume as defined above 
automatically satisfies the geometric conservation law (GCL), since from the 
Gauss theorem follows 

iidS = 0. (49) 



9CT 

In all numerical simulations shown later in this paper, it has been confirmed 
for all time steps and all elements that the above relation (49) was always 
satisfied on the discrete level up to machine precision. 

The final high order ALE one-step finite volume scheme takes the following 
form: 



1 1 



'n+l I r\n+l 



\Tr'\(^ 



TjdNi 



\dC2AG,^dxdT + 



c^\dc: 



Ph) d:sidt. 



(50) 

where the term Gjj-fijj contains the Arbitrary-Lagrangian-Eulerian numerical 
fiux function as well as the path-conservative jump term D, to resolve the 
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discontinuity of the predictor solution qh at the space-time sub-face dC^j. The 
surface integrals appearing in (50) are approximated using multidimensional 
Gaussian quadrature rules, see [126] for details. At the interface dC^j let us 
denote the local space-time predictor solution inside element Tj(t) by q^ and 
the element-local predictor solution of the neighbor element Tj{t) by q^, then 
a simple Rusanov-type scheme [56] is given by 



G^J = I (F(q+) + F(q,)) ■ fi., + ^ (/b(*) ■ n 



us U>max -I- 



^h-^h 



(51) 

where | Amax | is the maximum absolute value of the eigenvalues of the matrix 
A ■ n in space-time normal direction, which can be expressed in terms of 
the classical Eulerian system matrix A = dF/dQ + B and the normal mesh 
velocity V ■ n as 



A ■ n 



nl + n^ 



with 



n 



'OF 



B ■ n - (V ■ n) I 



(52) 



(53) 



As stated above, V • n is the local normal mesh velocity and I is the u x u 
identity matrix. It can be easily verified that 



/ 



V 




n, 



y^At 
-x^At 



\ 



hence V ■ n 



rit 



XyAjr Vx-^T I 



^nl + hi 

(54) 

A more sophisticated Osher-type scheme [110] has been introduced in the 
Eulerian framework for conservative and for non-conservative hyperbolic sys- 
tems in [54,55]. It has been extended to the Lagrangian framework in one 
space dimension in [56] and reads in the general multi-dimensional case with 
conservative and non-conservative terms as follows: 



G,, = 2(F(q+) + F(qj)-fi., + 2 



B(*l 



n 



An(*)|) ds\ (q^t-q^, 



(55) 
In (55) above, the usual definition of the matrix absolute value operator ap- 
plies, i.e. 

|A|=R|A|R-\ |A| =diag(|Ai|,|A2|,...,|A,|), (56) 

with the right eigenvector matrix R and its inverse R^^. According to [55,54] 
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the path integral appearing in (55) is approximated using Gaussian quadrature 
rules of sufficient accuracy. 



2.4 Numerical Convergence Studies 



In this section a numerical convergence study is performed for the compressible 
Baer-Nunziato model (1) in two space dimensions. This test problem has 
been proposed in [51] and has also been used in [55]. The test problem is 
similar to the one described in [76] and [8] . The exact solution of this smooth 
unsteady test problem is obtained in two steps: First, an exact stationary and 
rotationally symmetric solution of the governing PDE is sought and then the 
problem is made unsteady by superimposing a constant, uniform velocity field 
V using the principle of Galilean invariance of Newtonian mechanics. The exact 
solution is then simply given by the advection of the nontrivial initial condition 
with the superimposed constant velocity field v. The rotationally symmetric 
solution is found by writing the governing equations (1) in polar coordinates 
{r — 13) and by imposing angular symmetry (9/(9/3 = 0. What remains is an 
ODE system in the radial coordinate r that can be solved analytically, see 
[51]. 

In the following we denote with -uf the angular velocities and with -u^ the radial 
velocities. Since we are interested in a vortex-type solution, we furthermore 
suppose that u^ = 0. From the radial momentum equations we then obtain 
the following ODE system: 



2 

d ' ■ ° 1 / fl\ 

dr 



|:(0iPi) =P2|;0i + ^(mi) 0iPi, 

(02P2) = P2|;02 + I (^2) 02P2- 



If 01, Pi and P2 are know, e.g. by simply prescribing them, then (57) is just a 
simple algebraic equation system for the angular velocities wf . As in [51] we 
choose 

Pk=Pko (1 - ^e(^ """'/"Q] , (A; = 1,2) , (58) 



and 

3 2x/2^ 



. = 1 + -J=e-'-V2, (59) 



(60) 



hence the angular velocities of each phase result as 

u{ = ^^sjrD [pio (4v^Fi + QH^ - l2Gsl + SiJiS?) + 3p2oS? (4^ - H^ 

. r^/2 I — 

ul = JP2P20F2 , 

with the auxiliary variables 

_ 2r2 + r'^sl - 24 {r - Sk) {r + Sk) 

Hk = e 24 , Fk = e 4 , (A; = 1,2), 

and 

G' = e~"'/2^ D = Pi (2V2^ + 3G) . 

To this steady, rotationally symmetric solution of the compressible Baer- 
Nunziato equations we now add a constant uniform velocity field v = (m, v) 
to make the test problem unsteady, as already mentioned above. This can 
be done since Newtonian mechanics is Galilean invariant. With this manufac- 
tured analytical solution we can now calculate the convergence rates of the 
new class of high order Lagrangian one-step WENO finite volume schemes 
presented previously in this section. For the computational setup, we use the 
following parameters: 

7i = 1.4, 72 = 1.35, TTi = 7r2 = 0, u = v = 2, u = \ = 0, 

3 3 7 

Pi = 1, P2 = 2, pio = 1, P20 = -, Si = -, S2 = -. (61) 

The problem is solved with the Osher-type scheme (55) on a square domain 
Q = [—10; 10] X [—10; 10], using unstructured triangular meshes with four 
periodic boundary conditions, see Fig. 3, and setting the mesh velocity to 
V = u/ = Ui. The numerical convergence rates are shown for the solid volume 
fraction 0^ at time t = 2.0 in Table 1. One observes that the schemes reach 
their designed order of accuracy quite well. To our knowledge, this is the 
first time ever that a better than second order accurate Lagrangian WENO 
finite volume scheme is presented for non-conservative hyperbolic systems on 
unstructured triangular meshes with applications to the Baer-Nunziato model 
of compressible multi-phase fiows. 



3 Test Problems 



All the subsequent test problems are solved using the Osher-type method (55) 
and using the interface velocity, i.e. the solid velocity, as mesh velocity, hence 
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Table 1 

Numerical convergence results for the compressible Baer-Nunziato model using the 
third to sixth order version of the Arbitrary-Lagrangian-Eulerian one-step WENO 
finite volume schemes presented in this article. The error norms refer to the variable 
(j)s (solid volume fraction) at time t = 2.0. 



Ng 


eL2 


0{L2) 


Ng 


eL2 


0{L2) 




03 






04 




24 


2.6916E-02 


- 


24 


1.5993E-02 


- 


32 


1.0906E-02 


3.1 


32 


3.8281E-03 


5.0 


64 


1.9750E-03 


2.5 


64 


3.0900E-04 


3.6 


128 


2.5442E-04 


3.0 


128 


2.0855E-05 


3.9 




05 






06 




24 


1.4493E-02 


- 


24 


8.3869E-03 


- 


32 


3.8912E-03 


4.6 


32 


1.9504E-03 


5.1 


64 


2.5564E-04 


3.9 


64 


6.1843E-05 


5.0 


128 


8.7457E-06 


4.9 


96 


7.4509E-06 


5.2 



V = u/ = Ui. The CFL number in all test problems is set to CFL = 0.5. For 
all test problems we use a third order WENO scheme with reconstruction in 
characteristic variables [53] since componentwise reconstruction in conserva- 
tive variables led to significant spurious oscillations. In all cases shown below, 
friction and pressure relaxation are neglected, hence (A = z/ = 0). 



3.1 Riemann Problems 



The high order finite volume ALE schemes proposed in this article are subse- 
quently validated by applying them to ID Riemann problems that are solved 
in a 2D geometry on unstructured triangular meshes. The exact solution for 
these ID Riemann problems can be found in [5,120,39]. From the above men- 
tioned articles we have chosen a subset of four Riemann problems, whose initial 
conditions are listed in Table 2. Some of the test cases use the stiffened gas 
EOS, some of them consider just a mixture of two ideal gases. 

The initial two-dimensional computational domain is given by ^2(0) = [—0.5; 0.5] x 
[—0.05; 0.05], which is discretized using an initial characteristic mesh spacing 
oi h = 1/200, corresponding to an equivalent one-dimensional resolution of 
200 cells. The initial discontinuity is located at a; = and the final simulation 
times are listed in Table 2. In x-direction we use transmissive boundaries and 
in y-direction periodic boundary conditions are imposed. 

The numerical results are shown in Figs. 4-7 and are compared with the exact 
solution. On the top left of each figure a sketch of the mesh is depicted, while 
the other subfigures contain a one-dimensional cut through the reconstructed 
numerical solution w^ along the x-axis, evaluated at the final time on 200 
equidistant sample points. Due to the Lagrangian formulation of the method. 
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Fig. 3. Moving Lagrangian meshes used for the numerical convergence test at times 
t = (top), t = 1 (center) and t = 2 (bottom) with resolution 24 x 24 (left) and 
32 X 32 (right). 

the solid contact is resolved in a very sharp manner in all cases, which was 
actually the main aim in the design of a high order Lagrangian-type scheme 
for the compressible Baer-Nunziato model. Also for the other waves we can 
note in general a very good agreement between our numerical results and the 
exact reference solutions given in [5,120,39], apart for RP3, where a visible 
misfit between the exact solution and numerical solution is obtained for the 
gas pressure and the gas density. However, even a very diffusive central path- 
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conservative FORCE scheme in the Eulerian framework [51] had difficulties 
with this test problem and produced visible spurious oscillations and a wrong 
position of the solid contact in the gas phase. Further investigations on this 
problem are necessary, in particular, whether a different choice of the path is 
able to improve the situation. 



3.2 Cylindrical Explosion Problems 



Here we solve the compressible Baer-Nunziato equations in two space di- 
mensions on a circular domain Q{t) with initial radius R = 1.0. The initial 
condition is in all cases given by 

, -v., if |x| < 0.5, 

Q(x,0) = <{ - (62) 

else. 

The reference solution is obtained by solving an equivalent non-conservative 
one-dimensional PDE in radial direction with geometric reaction source terms, 
see [130] for the Euler equations and [131] for the Baer-Nunziato model for 
details. In our case here the reference solution has been obtained by using 
an Eulerian second order TVD scheme on a very fine ID mesh consisting of 
10,000 cells. The initial conditions for the three explosion problems solved here 
are taken from the Riemann problems solved previously, where the left state 
is used as inner state and the right state is used as the outer state in (62), 
respectively. In particular, the first explosion problem EPl uses the initial 
condition of RPl, EP2 corresponds to RP2 and EP3 to RP4, respectively. 
Also the parameters for the equation of state are chosen according to Table 
2. The initial mesh spacing is of characteristic size h = 1/250, leading to 

Table 2 

Initial states left (L) and right (R) for the Riemann problems solved in 2D and 3D 

with the Baer-Nunziato model. Values for 7j, ttj and the final time te are also given. 






Ps 


Us 


Ps Pg Ug Pg 


0s te 


RPl [39]: 


Is 


= 1.4, 


TTs = 0, 7g = 1.4, - 


/r3 = 


L 
R 


1.0 
2.0 


0.0 
0.0 


1.0 0.5 0.0 1.0 
2.0 1.5 0.0 2.0 


0.4 0.10 
0.8 


RP2 [39]: 


7s = 


= 3.0, 


vr, = 100, 7^ = 1.4, 


% = o 


L 
R 


800.0 
1000.0 


0.0 
0.0 


500.0 1.5 0.0 2.0 
600.0 1.0 0.0 1.0 
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0.3 


RP3 [39]: 


Is 


= 1.4, 


vTs = 0, 7g = 1.4, - 


Kg = Q 


L 
R 
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1.0 


0.9 
0.0 


2.5 1.0 0.0 1.0 
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0.2 


RP4 [120]: 


7s = 


3.0, 


vr^ = 3400, 7g = 1.35, 


% = o 


L 
R 


1900.0 
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10.0 2.0 0.0 3.0 
1000.0 1.0 0.0 1.0 


0.2 0.15 
0.9 
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Fig. 4. Results for Ricmann problem RPl of the seven-equation Baer-Nunziato 
model at time t = 0.1. 
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Fig. 5. Results for Ricmann problem RP2 of the seven-equation Baer-Nunziato 
model at time t = 0.1. 
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Fig. 6. Results for Ricmann problem RP3 of the seven-equation Baer-Nunziato 
model at time t = 0.1. 
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Fig. 7. Results for Riemann problem RP4 of the seven-equation Baer-Nunziato 
model at time t = 0.15. 

a total number of 431,224 triangular elements used to discretize Q{t). The 
numerical results are compared with the ID reference solution in Figures 8 
- 10. On the top left of each figure a 3D visualization of either the solid or 
the gas density is shown, in order to verify that the cylindrical symmetry 
is reasonably maintained on the unstructured triangular meshes used here. 
The other subfigures show a one-dimensional cut through the reconstructed 
numerical solution w/j on 250 equidistant sample points along the x-axis. In all 
cases an excellent agreement between numerical solution and reference solution 
is obtained. Note, in particular, the very sharp resolution of the solid contact 
due to the use of a Lagrangian framework. 



3.3 Two-Dimensional Riemann Problems 



In [83] Kurganov and Tadmor have collected a very nice set of numerical 
solutions for two-dimensional Riemann problems of the compressible Euler 
equations. Here, we propose an extension of these 2D Riemann problems to 
the compressible Baer-Nunziato model. The initial computational domain is 
Q[0) = [—0.5; 0.5] x [—0.5; 0.5] and the initial condition is given by four piece- 
wise constant states defined in each quadrant of the two-dimensional coordi- 
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Fig. 8. Results obtained for the first cylindrical explosion problem EPl at t = 0.15 
and comparison with the reference solution. 
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Fig. 9. Results obtained for the first cylindrical explosion problem EP2 at t = 0.15 
and comparison with the reference solution. 
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Fig. 10. Results obtained for the first cylindrical explosion problem EPS at t = 0.15 
and comparison with the reference solution. 
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Table 3 

Initial conditions for the two-dimensional Riemann problems. 





Configuration CI 


(7.= 


= 1.4, 


7, = 1-4 


,7r, = 


= TTg = 


= 0) 








Ps 


Us 


Vs 


Ps 


Pg 


Ug 


Vg 


Pg 


(ps 


Qi 


(x>0,y>0) 


2.0 


0.0 


0.0 


2.0 


1.5 


0.0 


0.0 


2.0 


0.8 


Q2 


(x<0,y>0) 


1.0 


0.0 


0.0 


1.0 


0.5 


0.0 


0.0 


1.0 


0.4 


Q3 


{x<0,y<0) 


2.0 


0.0 


0.0 


2.0 


1.5 


0.0 


0.0 


2.0 


0.8 


Q4 


{x>0,y<0) 


1.0 


0.0 


0.0 


1.0 


0.5 


0.0 


0.0 


1.0 


0.4 



Configuration C2 (7^ = 3.0, 7^ 

Ps Us Vs 



1.4, TT^ = 100, TTg = 0) 



Qi 

Q2 
Q3 
Q4 



(x > 0,y > 0) 
(x <0,y > 0) 
(x < 0,y < 0) 
{x > 0,y <0) 



Ps 



Ua 



Pg (h 
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800. 0.0 

1000. 0.0 

800. 0.0 
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0.0 
0.0 



600.0 
500.0 
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0.0 
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0.0 
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1.0 
2.0 



0.3 
0.4 
0.3 
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Q(x,0) 



fQi 
Q2 
Q3 

IQ4 



(63) 



nate system: 

if X > Ay > 0, 
if x<OAy>0, 
if X <OAy <0, 
if x>OAy<0. 

The initial conditions for the two configurations presented in this article are 
hsted in Table 3. The simulations are carried out with a third order one-step 
ALE WENO finite volume scheme using an unstructured triangular mesh 
composed of 90,080 elements with an initial characteristic mesh spacing of 
h = 1/200. The reference solution is computed with a high order Eulerian 
one-step scheme as presented in [131,55], using a very fine mesh composed 
of 2,277,668 triangles with characteristic mesh spacing h = 1/1000. Reflec- 
tive wall boundaries are applied on the four boundaries of the domain. The 
obtained results together with the Eulerian reference solution are depicted in 
Figures 12 - 13, where we can note a very good qualitative agreement of the 
Lagrangian solution with the Eulerian flne-grid reference solution. For the 
flrst test problem, the initial and the flnal mesh are depicted in Fig. 11. 
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Fig. 11. Mesh for configuration CI at times t = (left) and t = 0.15 (right). 
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Fig. 12. Results obtained with the third order Lagrangian WENO scheme for the 
2D Riemann problem CI at time t = 0.15 (left column). The reference solution 
computed with an Eulerian method on a very fine mesh is also shown (right column) . 
30 equidistant contour lines are shown for the solid density ps (top row), the gas 
density pg (middle row) and the solid volume fraction (ps (bottom row). 
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Fig. 13. Results obtained with the third order Lagrangian WENO scheme for the 
2D Riemann problem C2 at time t = 0.15 (left column). The reference solution 
computed with an Eulerian method on a very fine mesh is also shown (right column) . 
30 equidistant contour lines are shown for the solid density ps (top row), the gas 
density pg (middle row) and the solid volume fraction (ps (bottom row). 
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4 Application to Free Surface Flows with Moving Boundaries 



4-1 Reduced Three-Equation Baer-Nunziato Model 



The Baer-Nunziato model (1) considered in this paper can also be applied to 
complex non-hydrostatic free surface flows simulations, as proposed in [46,44]. 
There, a reduced three-equation Baer-Nunziato model has been considered, 
which can be derived from the complete seven-equation model (1) by in- 
troducing three assumptions: flrstly, all pressures are assumed to be relative 
pressures with respect to the atmospheric reference pressure, which is set to 
zero, i.e. po •= 0; secondly, the gas phase that surrounds the liquid is sup- 
posed to remain always at atmospheric reference conditions, hence obtaining 
P2 = Po = 0, which allows to neglect all evolution equations related to the 
gas phase; flnally, the pressure of the liquid phase is computed by the Tait 
equation of state [13], which reads 



Pi 




(64) 



where pi, po are the liquid density and the reference liquid density at standard 
conditions, respectively, ko is a constant that governs the compressibility of 
the fluid and 7 is a parameter commonly used to flt the equation of state with 
experimental data. Regarding the Tait equation (64), we should set the proper 
constant values for water, i.e. /cq = 3.2-10^ Pa, po = 1000 kg/rn? and 7 = 7, so 
that the typical speed of sound in water (c = 1500 m/ s) is preserved. However, 
since for most environmental free surface flows the maximum velocity is around 
10 m/ s, or even less, we would have to deal with extremely low Mach numbers 
of M = |u|/c ^ 1 that would cause very low time stepping and excessive 
numerical dissipation. In order to avoid this problem, we set artiflcially the 
Mach number to M = 0.3, as done in [46,44], and therefore we admit density 
fluctuations of the order of 10%. The advantage of using a weakly compressible 
model for free surface flows is the possibility to simulate environmental-type 
free surface flows as well as high speed industrial free surface flows, such as 
they appear in water jet cutting machines or fuel injection systems. There, 
speeds may reach up to |u| k. 1000 m/s and thus compressibility effects can 
no longer be neglected in the liquid. 

By using the above mentioned simpliflcations and introducing them into sys- 
tem (1), as outlined in [46,44], one can write the reduced Baer-Nunziato model 
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as 

^ (</)p) + V ■ (#u) = 0, 

d 

— (0pu) + V-(0(puu + cr)) = 0pg, (65) 

-0 + u-V0 = O, 

with the stress tensor cr = pi of the inviscid hquid phase, for which we have 
dropped the subscript 1 to ease notation. Mass and momentum equations are 
fully conservative in the system above, while the advection equation for the 
volume fraction is non-conservative. An evolution equation for total energy 
is not required, since we have used the assumption of a polytropic equation 
of state, which prescribes pressure as a function of density only. In (65) the 
state vector is Q = (0p, (ppu, 0) and g is the gravity vector acting along the 
vertical direction y, i.e. g = {0,g) with g = 9.81 m/s^. Further details on 
the derivation and on the assumptions of the reduced Baer-Nunziato model 
as well as a thorough validation against analytical solutions and experimental 
measurements can be found in [46,44]. The system (65) above involves a com- 
pressible inviscid phase, without viscous effects. However, viscosity plays an 
important role in turbulent flows with violent free surface motion and there- 
fore the algorithm can be improved and extended by adding the viscous terms 
to the stress tensor a of the momentum equation in (65). The stress tensor of 
a Newtonian fluid using the hypothesis of Stokes reads 

p + -/iV ■v^\- li (Vu + Vu'^) . (66) 

Here, /i denotes the dynamic viscosity, which for water is normally set to 10~^ 
N -s/rn?. The discretization of the above model (65), which now also contains 
the viscous effects, is done according to [45,48,72] and can be directly inserted 
into the high order one-step approach used in this paper. 



4-2 Sloshing in a Moving Tank 



We apply the reduced Baer-Nunziato model with viscosity (65) to a well 
known free surface flow problem in moving geometries, namely to the slosh- 
ing motion of a liquid in a partially fllled closed tank. Such a problem can 
not be described by the commonly used shallow water equations, since non- 
hydrostatic effects cannot be neglected. Furthermore, a numerical method is 
required that is able to deal with moving geometries, such as the present high 
order ALE flnite volume scheme. The flow is rather complex, characterized by 
the presence of high amplitude oscillations and, eventually, also wave breaking 
may occur. Since our algorithm is designed to be an Arbitrary Lagrangian- 
Eulerian (ALE) scheme, for the present problem we decide to move the mesh 
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according to the motion of the sloshing tank, which is prescribed on the bound- 
ary of the domain dQ{t). Inside the domain Q{t), the vertices are displaced 
smoothly according to the Laplace equation. In this case the evolution of the 
mesh and in particular the computation of the mesh velocity vector V(x, t) is 
governed by the following system of elliptic PDE: 

V2V(t) = 0, V(t) = Yoit) for X G dn{t), (67) 

which is solved in each time step and where the prescribed domain motion 
VD(t) is imposed as Dirichlet boundary condition. The solution is obtained 
using a classical PI finite element method (FEM), where the unknowns are 
located at the grid nodes and the solution of the discretized system (67) is 
computed by the conjugate gradient method. In this section the fluid motion 
and the mesh motion are governed by two different and independent laws, 
namely by the governing PDE system (65) and the Laplace equation (67), 
respectively. By solving (67) one obtains the velocity vector V for each grid 
vertex, which allows us to move the mesh nodes and to evolve the geometry, 
i.e. to compute all the other geometric quantities needed for the computation 
(normal vectors, element volumes, side lengths, barycenter positions, etc. ). 
From the solution of (67) the space-time control volumes needed for the local 
space-time Galerkin predictor and the final one-step finite volume scheme are 
known a priori. 

The motion of liquid sloshing is a classical and widely investigated problem, 
since it occurs in many real world conditions whenever a sway tank is present, 
as in cargo ships, in liquid tank carriages on rail roads or in propellant tanks 
of rockets and airplanes engines. Sloshing phenomena have gained recent at- 
tention in coastal and offshore engineering with the proliferation of liquefied 
natural gas (LNG) and oil carriers transporting liquids in partially filled tanks 
[139], in order to assure and guarantee the safety of the sea transport of those 
fuels. The liquid sloshing motion can become very complex, even including 
wave breaking phenomena that can create highly localized impact pressure on 
tank walls, hence causing structural damages and affecting the stability of the 
vehicle which carries the container. Sloshing phenomena have been widely and 
intensively investigated in the last decades, focusing on the development of an 
analytical study based on the potential flow theory, as done by Faltinsen [57] , 
who derived a linear analytical solution for liquid sloshing in a horizontally 
excited two-dimensional rectangular tank. Faltinsen and Timokha [59] extend 
the previous work to the nonlinear sloshing case and Hill [73] analyzed more 
in detail the waves' behavior by relaxing some of the simplifications intro- 
duced in the previous studies. Also laboratory measurements of wave height 
and hydrodynamic pressure have been collected and reported [135,104,105,4], 
which are very useful to validate both, theoretical solutions and numerical re- 
sults. The theoretical analyses are however not suitable to describe real fluid 
sloshing, where viscosity and turbulence occur, so that a lot of research has 
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been carried out in order to develop numerical models. In [57] Faltinsen de- 
veloped a boundary element method (BEM), while Nakayama and Washizu 
[103] analyzed the non-linear liquid sloshing in a two-dimensional rectangular 
tank under pitch excitation by using the finite element method (FEM). Wu 
et al. [137] were the first who conducted a series of three-dimensional demon- 
strations on liquid sloshing based on FEM, but they did not compare the 
results to any experimental data. Another numerical technique is given by the 
finite difference method (FDM) with the use of coordinate transformations, 
as proposed by Chen et al. [31], who adopted a curvilinear coordinate system 
to map the sloshing from the non-rectangular physical domain into a rectan- 
gular computational domain. One can also solve the Navier-Stokes equations 
for viscous liquid sloshing [29,28,30] or the Reynolds Averaged Navier-Stokes 
equations (RANS) as proposed by Armenio et al. [6], who observed that the 
RANS model provides much more accurate results than the classical shallow 
water equations (SWE). For liquid sloshing the vertical acceleration is indeed 
very important and can not be neglected. In literature one can also find La- 
grangian algorithms for sloshing phenomena, see the work of Okamoto et al. 
[104,105], who were the first to apply an Arbitrary Lagrangian-Eulerian finite 
element method to two- and three-dimensional liquid sloshing, and also mesh- 
less Lagrangian particle methods, such as smoothed particle hydrodynamics 
(SPH) algorithms are used, as proposed in [122]. Multi-fiuid sloshing has been 
studied with the Particle Finite Element Method (PFEM) in [79]. 

From the above mentioned references we know that the intensity and the 
impact of sloshing do generally depend on the amplitude and the frequency of 
the tank motion, on the physical properties of the fiuid and on the geometry 
of the tank, i.e. its shape. Furthermore, liquid sloshing occurs when a tank 
is partially filled with fiuid, so that a free surface is present. Very frequently 
high Reynolds numbers are achieved in such phenomena, therefore turbulence 
is usually present in this kind of flows. We improve indeed the algorithm by 
adding the simple zero-equation Smagorinsky turbulence model [123], which 
gives an expression for the eddy viscosity fir that is assumed to be proportional 
to the velocity gradients. For the two-dimensional case the Smagorinsky model 
reads 

I^T = P ■ {CsAf \Sl (68) 

where Cs is a coefficient which has to be set properly, A is the turbulence 
length scale and the term |S| is given by 




2S,,S,,, s^, = -U:^ + ^U (69) 



with the strain rate tensor Sij. The Smagorninsky constant is usually in the 
range [0.1; 0.24] and here we set Cs = 0.17. In this paper show numerical 
results for an idealized two-dimensional case, where the tank is moving with 
a purely horizontal sinusoidal velocity and the vertical velocity component is 
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Fig. 14. Sloshing tank. 



zero, I.e. 



yDit) = iVAsmicot),0f 



where u 



2tt 



(70) 

is the frequency of the oscillation and T is the period, while 
Va quantifies the amplitude of the sinusoidal velocity, which is related to the 
horizontal tank displacement 5xa by Va = —oo5xa- We always assume the 
fluid and the tank to be initially at rest, i.e. V(0) = 0. 

The simulations presented in this paper refer to the very recent work of Shao 
et al. [122], who present an improved SPH method for modeling liquid slosh- 
ing dynamics. Their main novelty regards the use of the Reynolds averaged 
turbulence model incorporated into the SPH method in order to properly 
describe the effects of turbulence. Furthermore, in [122] some new density 
and kernel gradient corrections have been introduced to achieve better ac- 
curacy and smoother pressure fields. The initial domain 0(0) refers to the 
one used for the laboratory measurements carried out by Faltinsen et al. 
[58] and is represented by a two-dimensional rectangular tank of dimensions 
r2(0) = [0; 1.73] X [0; 1.15], as sketched in Figure 14. We use an unstructured 
triangular mesh with a total number of A'"e = 44746 elements. The lateral sides 
and the bottom of the domain are modelled by classical no-slip wall boundary 
conditions for viscous flow, while the above boundary is set as transmissive 
outflow. The computational domain is characterized by a periodic excitation 
along the horizontal direction x, whose displacement 5x{t) is described as 
5x{t) = 6xa ■ cos(cjt), where we flx the parameters 6xa = 0.032 and T = 1.3 
with the initial water height hy^ = 0.6, as proposed in [122]. The corresponding 
velocity law for the tank can be easily computed by (70), with Va = —uj6xa- 
In this case we are dealing with a turbulent flow, since the Reynolds number 
Re = — is of the order of 10^. Regarding the parameters of the equation of 
state (64), we set ko = 8.5 • 10^ Pa, po = 1000 kg/rrv' and 7 = 1, as suggested 
in [46,44]. 

Figure 15 shows a comparison between experimental data and numerical re- 
sults: the periodic motion of the tank can be clearly identifled looking at the 
perturbation H of the free surface location with respect to the initial conflg- 



37 



uration. The numerical results have been collected at the same probe point 
Xp used for the experimental data, which is placed on the initial free sur- 
face level and is 0.05 m away from the left wall, hence its initial location is 
Xp(0) = (0.05,0.6). Since the probe point is attached to the tank, it moves 
together with the domain. Furthermore the value of H is evaluated at each 
time step as 

H = j (p{s)ds - K, with /if^ G [0,1.15] (71) 

where the volume fraction integral along the whole height h^, of the domain 
represents indeed the water column at the probe point which is shifting hori- 
zontally according to Eqn.(70) and h^ is the initial free surface elevation. As 
the tank motion begins, water moves towards the left, hence decreasing the 
pressure (t = 0.5), then the tank changes completely direction shifting towards 
the right side and therefore the free surface elevation increases its level up to 
0.7 m at time t = 1.2 s. The water keeps moving following the described peri- 
odic cycle of the tank boundaries, progressively increasing the wave amplitude 
in time. The numerical results are in good agreement with the experimental 
data, also with the particular wiggles observed in the time interval t G [5; 7]s, 
see Fig. 15. 

Another comparison of numerical results with experimental data is depicted 
in Figure 16 for two different sets of parameters: on the left we show the 
case T = 1.5s and hy^ = 0.6m, while on the right we use T = 1.875s and 
hw = 0.5m. The flow pattern of the sloshing for T = 1.5s and hy, = 0.6m is 
plotted in Figure 17 at nine representative time instants within one period. 
One can notice the strong oscillations of the free surface occurring in the flow 
while the tank is swinging, as well as the motion of the domain Q{t). 



5 Conclusions 



In this article the flrst better than second order Arbitrary-Lagrangian-Eulerian 
one-step WENO flnite volume scheme on unstructured triangular meshes has 
been proposed for the solution of hyperbolic systems with non-conservative 
products. The method has been applied to the full seven equation Baer- 
Nunziato model of compressible multiphase flows as well as to a reduced three- 
equation model for the simulation of weakly compressible free surface flows 
in moving geometries. High order of accuracy in space and time have been 
verifled by a numerical convergence study on a smooth unsteady test problem 
where an exact solution of the Baer-Nunziato model is available. The scheme 
has also been successfully applied to problems with shock waves and material 
interfaces. For all test problems the numerical results have been compared 
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Fig. 15. Sloshing motion: comparison between experimental data (solid line), second 
order (dotted line) and third order ALE WENO scheme (dash-dot line). T = 1.3 s 
and h^j = 0.6 m. 
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Fig. 16. Comparison between experimental data (solid line) and numerical results 
with second (dotted line) and third (dash-dot line) order of accuracy. Left: T = 1.5 
/i» = 0.5. Right: T = 1.875 /i^ = 0.5. 
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Fig. 17. Water motion for the case T = 1.5s and h^ 
3.4, 3.6, 3.8, 4.0, 4.2, 4.4 and 4.5 s. 



0.6m at times t = 3.0, 3.2, 



with exact or numerical reference solutions in order to validate the approach. 

In future research we plan to extend the present scheme to three space di- 
mensions in the more general framework of the new PnPm method proposed 
in [47], which unifies high order finite volume and discontinuous Galerkin fi- 
nite element methods in one general approach. Further work will also consist 
in a generalization to moving curved meshes as well as the introduction of 
multi-dimensional Riemann solvers, such as the ones used in [90,10,11]. Fi- 
nally, more research is necessary concerning a suitable conservative and high 
order accurate remapping / remeshing strategy to deal with significant mesh 
distorsions arising in the presence of strong shear waves. 
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